Pokaż/Ukryj kod
# Funkcja służąca do obliczenia prawdziwych wartości Y
Y_true <- function(x1, x2, epsilon){
return (0.5+x1+0.6*x2*x2+epsilon)
}Celem projektu jest zbadanie, w jaki sposób rozkład zmiennych objaśniających wpływa na wartość błędu średniokwadratowego ( \(MSE\) ) w modelu regresji liniowej.
Zagadnienie to jest istotne zarówno z teoretycznego, jak i praktycznego punktu widzenia - w szczególności w uczeniu maszynowym i analizie danych, gdzie poprawność założeń dotyczących rozkładu predyktorów wpływa na stabilność i trafność estymacji parametrów modelu.
Rozpatrujemy wpływ rozkładu zmiennych objaśniającyh na błąd średniokwadratowy (MSE) w następującym modelu regresji liniowej:
\[Y = 0.5 \ + \ X_1 \ + \ 0.6X_2^2 \ + \ \epsilon \]
gdzie \(X_1\) i \(X_2\) są zmiennymi losowymi pochodzącymi z rozkładów zależnych od rozpatrywanej sytuacji (patrz niżej), natomiast \(\epsilon\) jest szumem pochodzącym z rozkładu \(N(0,0.1)\) - jest tak, ponieważ chcemy się skupić wyłącznie na wpływie rozkładu predyktorów. MSE definiujemy następująco:
\[ MSE=\frac{1}{n} \sum_{i=1}^n(Y_i-\hat{Y_i})^2 \]
gdzie \(Y=(Y_1, Y_2, ..., Y_n)\) jest wektorem wartości obserwowanych, natomiast \(\hat{Y}=(Y_1, Y_2, ..., Y_n)\) jest wektorem wartości przewidywanych przez model, a \(n\) oznacza liczbę obserwacji w zbiorze danych. Dla całego projektu przyjmujemy, że \(n = 120\) , co pozwoli nam na bezpośrednie porównywanie uzyskiwanych wyników.
Poniżej przedstawiono wszystkie przypadki, w których generowane są zmienne losowe \(X_1\) i \(X_2\) . Opisano zarówno przypadki niezależne, jak i zależne, z uwzględnieniem możliwej zależności liniowej i nieliniowej tych zmiennych oraz wpływu zmiennej ukrytej na nie:
double_snorm
W tym przypadku zarówno \(X_1\) , jak i \(X_2\) pochodzą z rozkładu normalnego standaryzowanego: \(X_1\sim \mathcal{N}(0,1), \ X_2 \sim \mathcal{N}(0,1)\)
Ten przypadek jest naszym punktem odniesienia, ponieważ normalny rozkład standaryzowany jest symetryczny oraz ma umiarkowany rozrzut danych.
double_unif_narrow
Zarówno \(X_1\) , jak i \(X_2\) są losowane z wąskiego rozkładu jednostajnego, gdzie zakres danych jest ograniczony, a gęstości prawdopodobieństwa sa równomierne: \(X_1 \sim \mathrm{U}(-3,3), X_2 \sim \mathrm{U}(-3,3)\)
double_unif_wide
W tym przypadku zmienne są losowane z szerszego rozkładu jednostajnego: \(X_1\sim \mathrm{U}(-5,5), X_2 \sim \mathrm{U}(-5,5)\)
norm_exp
Zmienna \(X_1\) ma rozkład normalny z większym odchyleniem standardowym, a \(X_2\) pochodzi z rozkładu wykładniczego, co poszerza nasze badania o rozkład asymetryczny: \(X_1 \sim \mathcal{N}(0, 2.9), \quad X_2 \sim \mathrm{Exp}(\lambda = 1.2)\)
snorm_exp
Tutaj \(X_1\) ma rozkład normalny standaryzowany, a \(X_2\) - wykładniczy. Przypadek ten w zestawieniu z przypadkiem 4. pozwala zobaczyć wpływ wariancji rozkładu normalnego na MSE: \(X_1 \sim \mathcal{N}(0,1), \quad X_2 \sim \mathrm{Exp}(\lambda = 1.2)\)
bimodal_snorm
W tym przypadku \(X_1\) ma rozkład bimodalny, a \(X_2\) - normalny standaryzowany: \(X_1 \sim 0.43\cdot \mathcal{N}(-3,0.8) + 0.57\cdot \mathcal{N}(3,0.8), \quad X_2 \sim \mathcal{N}(0,1)\)
Dzięki temu zobaczymy, jak zmienia się MSE w przypadku gdy mamy dwie mody w jednym z rozkładów.
exp_norm
Zmienna \(X_1\) pochodzi z rozkładu wykładniczego, a \(X_2\) z normalnego: \(X_1 \sim \mathrm{Exp}(\lambda = 1.2), \quad X_2 \sim \mathcal{N}(0, 2.9)\)
Tutaj sprawdzamy sytuację w pewnym sensie odwrotną do przypadku 4., czyli co się dzieje, gdy to zmienna \(X_1\) ma duży ogon rozkładu.
t_snorm
Zmienna \(X_1\) pochodzi z rozkładu t-Studenta z 3 stopniami swobody, a \(X_2\) z rozkładu normalnego standaryzowanego: \(X_1 \sim t(df=3), \quad X_2 \sim \mathcal{N}(0,1)\)
t_norm
Zmienna \(X_1\) pochodzi z rozkładu t-Studenta, a \(X_2\) z normalnego rozkładu o większej wariancji: \(X_1 \sim t(df=3), \quad X_2 \sim \mathcal{N}(0,2.9)\)
t_exp
Zmienna \(X_1\) pochodzi z rozkładu t-Studenta, natomiast \(X_2\) z rozkładu wykładniczego: \(X_1 \sim t(df=3), \quad X_2 \sim \mathrm{Exp}(\lambda = 1.2)\)
W przypadkach 10.-12. analizujemy wpływ grubych ogonów i wartości odstających na MSE. Te trzy sytuacje możemy porównywać z dowolnym z przypadków 1-9, gdzie zmienna objaśniająca \(X_2\) ma dokładnie taki sam rozkład jak w przypadkach 10-12.
Korelacja liniowa
linear_corr_snorm Obie zmienne pochodzą z rozkładu normalnego standaryzowanego, ale wprowadzona jest między nimi korelacja \(\rho\) : \((X_1, X_2) \sim \mathcal{N}\Big(\mu =\begin{pmatrix}0\\0\end{pmatrix},\Sigma =\begin{pmatrix}1 & \rho \\ \rho & 1\end{pmatrix}\Big)\)
Taki przypadek pozwala badać wpływ liniowej zależności między predyktorami
Zależność nieliniowa nonlinear_corr Zmienna \(X_2\) pochodzi z rozkładu normalnego standaryzowanego. Zmienna \(X_1\) jest nieliniową funkcją zmiennej \(X_2\) z dodatkiem szumu \(\eta\) : \(X_2 \sim \mathcal{N}(0,1), \quad X_1 = a \cdot e^{b X_2} + \eta, \quad \eta \sim \mathcal{N}(0,0.1)\)
Pozwala to na modelowanie przypadków, w których zależność między predyktorami jest wyraźnie nieliniowa
Zmienne zależne od zmiennej ukrytej Z hidden_impact
Tutaj zarówno \(X_1\) , jak i \(X_2\) zależą od wspólnej zmiennej ukrytej \(Z\) : \(X_1 = 2.7 Z + 12, \quad X_2 = 1.4 Z + 15, \quad Z \sim \mathcal{N}(0,1)\)
Taki przypadek umożliwia zbadanie wpływu zmiennej ukrytej na obserwowane zmienne objaśniające
Razem zbadamy 17 przypadków. Dla każdego z nich generować będziemy po 200 ramek danych, co da nam wystarczającą próbkę do porównywania wyników. Dla każdej z utworzonych ramek danych będziemy dopasowywać model regresji liniowej, zapiszemy uzyskane wyniki i porównamy je. Tym samym wyciągniemy najciekawsze informacje z wpływu różnych rozkładów predyktorów na wartość błędu średniokwadratowego tworzonej regresji.
Wprowadzamy funkcję Y_true , która oblicza wartości wektora Y na podstawie wygenerowanych danych. Ponadto dodajemy epsilon oraz eta, czyli szum, aby nadać prawdziwym wartościom \(Y\) losowości. W ten sposób tworzymy “prawdziwe, rzeczywiste” dane.
Do generowania danych do każdego z wymienionych we wstępie przypadków rozkładu danych zmiennych objaśniających wykorzystujemy funkcję generate_data . Uwzględnia ona wszystkie założenia rozpatrywanych przez nas w tym projekcie scenariuszy. Generuje odpowiednie dane w zależności od wprowadzonego argumentu case (co tłumaczymy jako ‘przypadek’). Następnie wygenerowane dane są zapisywane do ramki danych, która zawiera 4 kolumny:
X1 - wartości pierwszej zmiennej objaśniającej
X2 - wartości drugiej zmiennej objaśniającej
Y - obliczone na podstawie wzoru (numer wzoru) wartości Y
Eps - biały szum
Funkcja zwraca tę ramkę danych oraz meta dane generatora rozkładu, czyli listę parametrów użytych podczas generowania danych. Parametry te można poznać poniżej. A zatem wynik funkcji generate_data to lista posiadająca dwa elementy:
case# Funkcja służąca do obliczenia prawdziwych wartości Y
Y_true <- function(x1, x2, epsilon){
return (0.5+x1+0.6*x2*x2+epsilon)
}generate_data <- function(case,
n = 120, # liczba generowanych obserwacji
rho = NULL, # współczynnik korelacji (przypadek zależności liniowej)
lambda = 2, # lambda dla rozkładu Poisson'a
exp_rate = 1.2, # intensywność dla rozkładu wykładniczego
a = NULL, # wartość pierwszego współczynnika dla zależności nieliniowej
b = NULL # wartość drugiego współczynnika dla zależności nieliniowej
){
# definicja szumu
epsilon <- rnorm(120, sd = 0.1)
# Globalna definicja szumu dla przypadku zmiennej zależnej
eta <- rnorm(120, sd = 0.1)
# sprawdzamy, czy został podany obowiązkowy parametr dla danego przypadku generowania danych
needed <- function(param, name){
if(is.null(param)){
stop(sprintf("Przypadek '%s' wymaga argumentu '%s'", case, name))
}
}
# przypadek, gdy oba predyktory mają rozkład normalny standaryzowany
if(case=="double_snorm"){
x1 <- rnorm(n)
x2 <- rnorm(n)
# Zmienne X_1 i X_2 pochodzą z rozkładu jednostajnego o parametrach (-3,3)
}else if(case == "double_unif_narrow"){
x1 <- runif(n, min = -3, max = 3)
x2 <- runif(n, min = -3, max = 3)
# Zmienne X_1 i X_2 pochodzą z rozkładu jednostajnego o parametrach (-5,5)
}else if(case == "double_unif_wide"){
x1 <- runif(n, min = -5, max = 5)
x2 <- runif(n, min = -5, max = 5)
# Zmienna X_1 pochodzi z rozkładu normalnego z większym odchyleniem standardowym,
# natomiast zmienna X_2 pochodzi z rozkładu wykładniczego z lambdą = 1.2 (parametr intensywności)
}else if(case == "norm_exp"){
needed(exp_rate,"exp_rate")
x1 <- rnorm(n, sd = 2.9)
x2 <- rexp(n,rate = exp_rate)
# Zmienna X_1 pochodzi z rozkładu normalnego standaryzowanego, a zmienna X_2
# z rozkładu wykładniczego o lambdzie równej 1.2
}else if(case == "snorm_exp"){
needed(exp_rate,"exp_rate")
x1 <- rnorm(n)
x2 <- rexp(n, rate = exp_rate)
# Zmienna X_1 pochodzi z rozkładu bimodalnego, a zmienna X_2 z normalnego
# standaryzowanego
}else if(case == "bimodal_snorm"){
bimod_mean1 <- -3
bimod_mean2 <- 3
bimod_sd <- 0.8
p <- 0.43
u <- runif(n)
x1 <- ifelse(u<p, rnorm(n, bimod_mean1, bimod_sd),
rnorm(n,bimod_mean2,bimod_sd))
x2 <- rnorm(n)
#Przypadek, gdy zmienna X_1 pochodzi z rozkładu wykładniczego,
# a zmienna X_2 z rozkładu normalnego
}else if(case == "exp_norm"){
needed(exp_rate,"exp_rate")
x1 <- rexp(n, rate = exp_rate)
x2 <- rnorm(n, sd = 2.9)
}
#X_1 pochodzi z rozkładu t-Studenta o trzech stopniach swobody,
# a X_2 z rozkładu normalnego standaryzowanego
else if(case == "t_snorm"){
x1 <- rt(n,df = 3)
x2 <- rnorm(n)
#X_1 pochodzi z rozkładu t-Studenta o trzech stopniach swobody,
# a X_2 z rozkładu normalnego
}else if(case == "t_norm"){
x1 <- rt(n, df = 3)
x2 <- rnorm(n,sd = 2.9)
#X_1 pochodzi z rozkładu t-Studenta o trzech stopniach swobody,
# a X_2 z rozkładu wykładniczego
}else if(case == "t_exp"){
needed(exp_rate,"exp_rate")
x1 <- rt(n, df = 3)
x2 <- rexp(n, rate = exp_rate)
# Zmienne losowe X_1 i X_2 są liniowo skorelowane i pochodzą obie
# z rozkładu normalnego standaryzowanego
}else if(case == "linear_corr_snorm"){
needed(rho,"rho")
sigma <- rbind(c(1,rho), c(rho,1)) #macierz kowariancji
mu <- c(0,0) #średnie rozkładów
# Zmienna X_1 jest nieliniowo zależna od zmiennej X_2
# a zmienna X_2 pochodzi z rozkładu normalnego standaryzowanego
}else if(case == "nonlinear_corr"){
needed(a,"a")
needed(b,"b")
x2 <- rnorm(n)
x1 <- a*exp(b*x2)+eta
}else if(case == "hidden_impact"){
z <- rnorm(n)
x1 <- 2.7*z + 12
x2 <- 1.4*z + 15
}
meta <- list(case = case, n = n,
params = list(rho = rho, lambda = lambda, exp_rate = exp_rate))
if(case == "linear_corr_snorm"){
df <- as.data.frame(mvrnorm(n=n, mu=mu, Sigma=sigma))
df <- df %>%
mutate(Y=Y_true(V1, V2, epsilon),
Eps = epsilon)
colnames(df) <- c("X1","X2","Y", "Eps")
return(list(df=df,meta=meta))
}else if(case == "hidden_impact"){
df <- data.frame(X1 = x1, X2 = x2, Y = 0.5+x1+0.6*x2*x2+z+epsilon, Eps=epsilon)
colnames(df) <- c("X1","X2","Y", "Eps")
return(list(df = df, meta = meta, z = z))
}
y <- Y_true(x1,x2,epsilon)
df <- data.frame(X1 = x1, X2 = x2, Y = y, Eps = epsilon)
colnames(df) <- c("X1","X2","Y", "Eps")
return(list(df=df,meta=meta))
}Sprawdźmy działanie naszych funkcji na poszczególnych przypadkach:
Uwaga! Funkcja ggpairs wykonuje tzw. kernel density estimation podczas rysowania funkcji gęstości na podstawie realizacji prób pochodzących z przyjętych teoretycznych rozkładów. Biorąc pod uwagę również to, że nie generujemy bardzo dużej próbki danych dla naszych zmiennych, otrzymujemy przybliżone wykresy gęstości dla rozkładów prawdopodobieństwa ze wstępu na przekątnych widzianych poniżej obrazków.
Poniższa tabela zawiera empiryczne i teoretyczne informacje dotyczące wygenerowanych ramek danych. Pozwala to na porównanie poprawności otrzymywanych rozkładów.
| Rozklad | mean_X1 | expected_mean_x1 | var_X1_mean | expected_var_x1 | mean_X2 | expected_mean_x2 | var_X2_mean | expected_var_x2 |
|---|---|---|---|---|---|---|---|---|
| double_snorm | 0.0116 | 0.0000 | 0.9948 | 1.0000 | 0.0063 | 0.0000 | 0.9918 | 1.000 |
| double_unif_narrow | 0.0064 | 0.0000 | 3.0261 | 3.0000 | -0.0045 | 0.0000 | 2.9869 | 3.000 |
| double_unif_wide | 0.0021 | 0.0000 | 8.2778 | 8.3330 | 0.0172 | 0.0000 | 8.4556 | 8.333 |
| norm_exp | -0.0253 | 0.0000 | 8.7116 | 8.4100 | 0.8311 | 0.8333 | 0.6993 | 0.694 |
| snorm_exp | 0.0009 | 0.0000 | 1.0516 | 1.0000 | 0.8262 | 0.8333 | 0.6921 | 0.694 |
| bimodal_snorm | 0.4722 | 0.4200 | 9.5479 | 9.4636 | -0.0171 | 0.0000 | 1.0165 | 1.000 |
| exp_norm | 0.8243 | 0.8333 | 0.6496 | 0.6940 | -0.0072 | 0.0000 | 8.4423 | 8.410 |
| t_snorm | 0.0332 | 0.0000 | 2.8118 | 3.0000 | 0.0272 | 0.0000 | 1.0145 | 1.000 |
| t_norm | -0.0044 | 0.0000 | 2.7619 | 3.0000 | 0.0113 | 0.0000 | 8.5536 | 8.410 |
| t_exp | -0.0480 | 0.0000 | 3.2950 | 3.0000 | 0.8366 | 0.8333 | 0.7186 | 0.694 |
| linear_corr_snorm | 0.0081 | 0.0000 | 1.0332 | 1.0000 | 0.0055 | 0.0000 | 1.0400 | 1.000 |
| nonlinear_corr | 3.2725 | 3.0802 | 116.2893 | 80.5394 | 0.0208 | 0.0000 | 0.9905 | 1.000 |
| hidden_impact | 11.9319 | 0.0000 | 7.3110 | 7.3000 | 14.9647 | 0.0000 | 1.9656 | 1.970 |
Każdy wiersz powyższej tabeli odpowiada jednemu z przypadków symulacyjnych, a każda kolumna przedstawia empiryczne i teoretyczne parametry zmiennych \(X_1\) i \(X_2\) obliczone na podstawie 30 powtórzeń generowania ramek danych z tymi zmiennymi oraz na podstawie teoretycznych obliczeń odpowiednio.
Rozklad – nazwa przypadku symulacyjnego, np. t_norm, pois_snorm, bimodal_snorm itd.
mean_X1 – średnia wartość empiryczna zmiennej \(X_1\) , obliczona jako średnia ze średnich z 30 wygenerowanych prób
expected_mean_x1 – teoretyczna wartość oczekiwana zmiennej \(X_1\)
var_X1_mean – średnia z oszacowanych wariancji \(X_1\) z 30 ramek danych
expected_var_x1 – teoretyczna wariancja zmiennej \(X_1\)
mean_X2 – średnia empiryczna zmiennej \(X_2\) (średnia ze średnich z 30 wygenerowanych ramek danych)
expected_mean_x2 – teoretyczna wartość oczekiwana zmiennej \(X_2\)
var_X2_mean – średnia z oszacowanych wariancji \(X_2\) z 30 ramek danych
expected_var_x2 – teoretyczna wariancja zmiennej \(X_2\)
Przypominamy raz jeszcze, że ta sekcja dotyczyła wyłącznie testowania funkcji do generowania danych. Uzyskane w tej sekcji wyniki dają nam praktyczny wgląd na sposób generowania danych przez dedykowaną funkcję. Widzimy, że rozkłady zgadzają się z założeniami, a tworzone ramki danych są miarodajne. Możemy zatem przejść do dalszych badań.
lm()?Model lm() w R implementuje klasyczną metodę najmniejszych kwadratów (MNK), która znajduje takie współczynniki \(\beta\), aby minimalizować sumę kwadratów błędów między wartościami obserwowanymi (Y) a wartościami przewidywanymi przez model \((\hat{Y})\).
Oznacza to, że lm() szuka prostych (hiperpłaszczyzn) najlepiej dopasowanych do danych w sensie średniokwadratowym (czyli właśnie tego, co mierzy MSE).
Funkcja prawdziwa postaci \(Y=0.5+X_1+0.6*X_2^2 +\epsilon\) jest modelem liniowym względem parametrów, ponieważ zależność od \(\beta_0\), \(\beta_1\), \(\beta_2\) jest liniowa (nawet jeśli występuje kwadrat zmiennej \(X_2^2\)). Dlatego możemy używać lm() nawet mimo obecności nieliniowości w zmiennych.
Wybraliśmy trzy różne modele aby ocenić, jak różny stopień dopasowania modelu do prawdziwej zależności wpływa na wynik metryki MSE.
Każdy z modeli reprezentuje inny poziom trafności założeń:
Model Ideal: lm(Y ~ X1 + I(X2^2)) odzwierciedla faktyczną postać funkcji generującej dane (najlepsze możliwe dopasowanie)
Model Simple: lm(Y ~ X1 + X2) jest celowo uproszczony, aby pokazać, jak wzrasta błąd, gdy pominiemy ważny nieliniowy element - służy jako punkt odniesienia.
Model Mixed: lm(Y ~ X1 + X2 + I(X1*X2)) model bardziej złożony niż klasyczny model liniowy
I() mówi funkcji lm() lub innej modelującej funkcji, że to, co jest w środku, ma być zinterpretowane dosłownie jako wyrażenie matematyczne, a nie jako specjalna składnia formuły R.
library(dplyr)
library(ggplot2)
library(MASS)
library(broom)
library(gridExtra)
set.seed(123)
cases <- c(
"double_snorm",
"double_unif_narrow",
"double_unif_wide",
"norm_exp",
"snorm_exp",
"bimodal_snorm",
"exp_norm",
"t_snorm",
"t_norm",
"t_exp",
"linear_corr_snormrho0.2",
"linear_corr_snorm_rho0.5",
"linear_corr_snorm_rho0.9",
"nonlinear_corr_a0.5b0.2",
"nonlinear_corr_a1b0.5",
"nonlinear_corr_a1.5b1",
"hidden_impact"
)
n_iter <- 200
n_cases <- length(cases)
n_obs <- 120# Zapis przykładowych ramek danych (po jednej dla przykładu)
if (!dir.exists("data_samples")) dir.create("data_samples")
for (case_name in cases) {
if (grepl("linear_corr_snorm", case_name)) {
rho_val <- as.numeric(sub(".*rho", "", case_name))
df_list <- generate_data(case = "linear_corr_snorm", rho = rho_val)
} else if (grepl("nonlinear_corr", case_name)) {
ab_vals <- regmatches(case_name, gregexpr("[0-9\\.]+", case_name))[[1]]
a_val <- as.numeric(ab_vals[1])
b_val <- as.numeric(ab_vals[2])
df_list <- generate_data(case = "nonlinear_corr", a = a_val, b = b_val)
} else {
df_list <- generate_data(case = case_name)
}
df <- df_list$df
write.csv(df, paste0("data_samples/", case_name, ".csv"), row.names = FALSE)
}Każdy zbiór danych jest dzielony na:
80% train - do dopasowania modelu (lm())
20% test - do oceny jakości dopasowania (liczymy MSE)
Jeśli policzylibyśmy błąd MSE na tych samych danych, na których model był trenowany, mielibyśmy zaniżony błąd (overfitting).
Podział pozwala realnie ocenić zdolność modelu do generalizacji, czyli jak radzi sobie z nowymi danymi.
W taki sposób podzielone zostały wszystkie ramki danych (200) dla każdego przypadku.
MSE wyraża rozproszenie danych od linii regresji.
Mniejsze MSE oznacza, że model jest lepiej dopasowany, błędy predykcji są mniejsze
Większe MSE oznacza, że model gorzej odwzorowuje zależności w danych.
Należy być jednak ostrożnym, ponieważ MSE jest miarą bezwzględną - w zależności od sytuacji znaczenie wielkości może się zmieniać. Przy analizie konkretnych przypadków będziemy zwracać uwagę również na wariancję zmiennych, która wpływa na wartość błędu średniokwadratowego.
MSE nie jest bezpośrednio porównywalne między przypadkami, jeśli skale \(X_{1}\) i \(X_2\) różnią się.
Przykładowo dane o dużym rozrzucie np. \(X\sim N(0,10)\) mogą mieć naturalnie większe błędy prognozy niż dane z \(X \sim N(0,1)\), nawet jeśli dopasowanie względne jest podobne. Dlatego przy analizie należy wziąć pod uwagę wariancję \(X\).
# dopasowanie modeli i liczenie MSE
results_list <- list()
idx <- 1
for (case_id in 1:n_cases) {
case_name <- cases[case_id]
for (i in 1:n_iter) {
# generowanie danych
if (grepl("linear_corr_snorm", case_name)) {
# wyciągamy rho z nazwy, np. z 'rho0.5'
rho_val <- as.numeric(sub(".*rho", "", case_name))
df_list <- generate_data(case = "linear_corr_snorm", rho = rho_val)
} else if (grepl("nonlinear_corr", case_name)) {
# wyciągamy a i b z nazwy, np. z 'a1.5b1'
ab_vals <- regmatches(case_name, gregexpr("[0-9\\.]+", case_name))[[1]]
a_val <- as.numeric(ab_vals[1])
b_val <- as.numeric(ab_vals[2])
df_list <- generate_data(case = "nonlinear_corr", a = a_val, b = b_val)
} else {
df_list <- generate_data(case = case_name)
}
df <- df_list$df
# Split 80/20
train_idx <- sample(1:nrow(df), size = 0.8 * nrow(df))
train <- df[train_idx, ]
test <- df[-train_idx, ]
m_simple <- lm(Y ~ X1 + X2, data = train)
m_ideal <- lm(Y ~ X1 + I(X2^2), data = train)
m_mixed <- lm(Y ~ X1 + X2 + I(X1 * X2), data = train)
# wyliczanie MSE i zapis współczynników
models <- list(simple = m_simple, ideal = m_ideal, mixed = m_mixed)
for (m_name in names(models)) {
m <- models[[m_name]]
preds <- predict(m, newdata=test)
mse <- mean((test$Y - preds)^2)
coefs <- coef(m)
row <- tibble(
case = case_name,
iter = i,
model = m_name,
mse = mse,
intercept = ifelse("(Intercept)" %in% names(coefs), coefs["(Intercept)"], NA_real_),
b_X1 = ifelse("X1" %in% names(coefs), coefs["X1"], NA_real_),
b_X2 = ifelse("X2" %in% names(coefs), coefs["X2"], NA_real_),
b_X2sq = ifelse("I(X2^2)" %in% names(coefs), coefs["I(X2^2)"], NA_real_),
b_X1X2 = ifelse("I(X1 * X2)" %in% names(coefs), coefs["I(X1 * X2)"], NA_real_)
)
results_list[[idx]] <- row
idx <- idx + 1
}
}
}
results_all <- bind_rows(results_list)
# Zapis wyników
write.csv(results_all, "wyniki_parametry_MSE.csv", row.names = FALSE)MSE jest zmienną losową, ponieważ zależy od losowo wygenerowanych danych. Powtarzamy eksperyment 200 razy dla każdego przypadku, co pozwala nam policzyć
średnie - przeciętny błąd
mediany - błąd “typowy”, bardziej odporny na wartości odstające
kwartyle Q25 i Q75, które pokazują rozrzut błędów
odchylenia standardowe - mówi o zmienności wyników między iteracjami
i ocenić stabilność modelu - czy między iteracjami MSE się waha, czy jest stabilne. Sprawdzamy w ten sposób w jakich rozkładach błędy są największe i jak duże znaczenie w dopasowaniu się do danych miał dobór wzoru modelu.
# Statystyki MSE
results_summary <- results_all %>%
group_by(case, model) %>%
summarise(
min_MSE=min(mse),
q25 = quantile(mse, 0.25),
mean_MSE = mean(mse),
median_MSE = median(mse),
q75 = quantile(mse, 0.75),
max_MSE=max(mse),
sd_MSE = sd(mse),
.groups = "drop"
)
write.csv(results_summary, "summary_MSE.csv", row.names = FALSE)
results_summary %>%
mutate(across(where(is.numeric), ~round(., 4))) %>%
kable(caption = "Statystyki MSE dla poszczególnych przypadków i modeli") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
column_spec(1, bold = TRUE) %>%
collapse_rows(columns = 1, valign = "middle")| case | model | min_MSE | q25 | mean_MSE | median_MSE | q75 | max_MSE | sd_MSE |
|---|---|---|---|---|---|---|---|---|
| bimodal_snorm | ideal | 0.0042 | 0.0082 | 0.0102 | 0.0096 | 0.0120 | 0.0189 | 0.0029 |
| mixed | 0.1354 | 0.4141 | 0.7973 | 0.6696 | 0.9658 | 3.7760 | 0.5559 | |
| simple | 0.1473 | 0.4101 | 0.7559 | 0.5939 | 0.9418 | 4.1909 | 0.5317 | |
| double_snorm | ideal | 0.0035 | 0.0083 | 0.0104 | 0.0102 | 0.0126 | 0.0191 | 0.0030 |
| mixed | 0.1561 | 0.4176 | 0.8112 | 0.6740 | 1.0470 | 2.8149 | 0.5103 | |
| simple | 0.1596 | 0.4133 | 0.7927 | 0.6619 | 1.0557 | 3.1058 | 0.5136 | |
| double_unif_narrow | ideal | 0.0034 | 0.0085 | 0.0107 | 0.0105 | 0.0128 | 0.0188 | 0.0029 |
| mixed | 1.2042 | 2.3274 | 2.7298 | 2.6784 | 3.0705 | 4.6891 | 0.6333 | |
| simple | 1.2037 | 2.2756 | 2.6996 | 2.6607 | 3.0154 | 4.5204 | 0.6129 | |
| double_unif_wide | ideal | 0.0028 | 0.0083 | 0.0105 | 0.0103 | 0.0125 | 0.0221 | 0.0032 |
| mixed | 9.2860 | 18.1913 | 21.7211 | 21.5384 | 25.1638 | 39.0761 | 5.4007 | |
| simple | 9.8518 | 18.0734 | 21.2760 | 21.6044 | 24.2744 | 37.9374 | 4.9723 | |
| exp_norm | ideal | 0.0047 | 0.0080 | 0.0102 | 0.0101 | 0.0122 | 0.0196 | 0.0029 |
| mixed | 9.7640 | 26.9976 | 55.1848 | 45.8961 | 62.8221 | 315.2500 | 46.9501 | |
| simple | 9.1819 | 25.8804 | 52.6158 | 43.2131 | 62.5845 | 314.3302 | 42.3956 | |
| hidden_impact | ideal | 0.0041 | 0.0085 | 0.0105 | 0.0101 | 0.0121 | 0.0214 | 0.0031 |
| mixed | 0.0041 | 0.0085 | 0.0105 | 0.0101 | 0.0121 | 0.0214 | 0.0031 | |
| simple | 0.4090 | 1.5255 | 3.1123 | 2.4248 | 3.8655 | 14.7734 | 2.2301 | |
| linear_corr_snorm_rho0.5 | ideal | 0.0043 | 0.0080 | 0.0102 | 0.0098 | 0.0122 | 0.0201 | 0.0031 |
| mixed | 0.1162 | 0.2725 | 0.4740 | 0.3749 | 0.5575 | 2.0048 | 0.3348 | |
| simple | 0.1394 | 0.3853 | 0.7723 | 0.5959 | 0.9238 | 3.7862 | 0.5961 | |
| linear_corr_snorm_rho0.9 | ideal | 0.0044 | 0.0083 | 0.0105 | 0.0105 | 0.0128 | 0.0192 | 0.0029 |
| mixed | 0.0215 | 0.0587 | 0.0904 | 0.0815 | 0.1144 | 0.2518 | 0.0448 | |
| simple | 0.1730 | 0.3888 | 0.7513 | 0.6216 | 0.9944 | 2.6078 | 0.4840 | |
| linear_corr_snormrho0.2 | ideal | 0.0043 | 0.0084 | 0.0104 | 0.0103 | 0.0122 | 0.0187 | 0.0027 |
| mixed | 0.1561 | 0.4055 | 0.7606 | 0.5901 | 0.9027 | 3.1477 | 0.5390 | |
| simple | 0.1749 | 0.4532 | 0.7926 | 0.6427 | 0.9488 | 3.3684 | 0.5357 | |
| nonlinear_corr_a0.5b0.2 | ideal | 0.0036 | 0.0080 | 0.0102 | 0.0098 | 0.0121 | 0.0198 | 0.0029 |
| mixed | 0.0552 | 0.1596 | 0.2614 | 0.2192 | 0.3074 | 1.4891 | 0.1718 | |
| simple | 0.1673 | 0.4194 | 0.7529 | 0.6097 | 0.8760 | 4.5577 | 0.5439 | |
| nonlinear_corr_a1.5b1 | ideal | 0.0042 | 0.0081 | 0.0102 | 0.0098 | 0.0121 | 0.0228 | 0.0031 |
| mixed | 0.0202 | 0.0415 | 0.0751 | 0.0563 | 0.0749 | 0.8991 | 0.0802 | |
| simple | 0.0200 | 0.0911 | 0.2847 | 0.1427 | 0.2493 | 4.0425 | 0.5114 | |
| nonlinear_corr_a1b0.5 | ideal | 0.0047 | 0.0086 | 0.0104 | 0.0097 | 0.0121 | 0.0187 | 0.0028 |
| mixed | 0.0154 | 0.0553 | 0.1458 | 0.0937 | 0.1452 | 2.1437 | 0.2036 | |
| simple | 0.0601 | 0.1383 | 0.2110 | 0.1764 | 0.2370 | 1.9686 | 0.1618 | |
| norm_exp | ideal | 0.0037 | 0.0078 | 0.0102 | 0.0100 | 0.0121 | 0.0272 | 0.0034 |
| mixed | 0.0537 | 0.1902 | 0.8671 | 0.3668 | 0.6923 | 19.8119 | 2.0971 | |
| simple | 0.0522 | 0.1951 | 0.7535 | 0.3549 | 0.6682 | 15.8515 | 1.6865 | |
| snorm_exp | ideal | 0.0045 | 0.0083 | 0.0100 | 0.0098 | 0.0114 | 0.0188 | 0.0026 |
| mixed | 0.0441 | 0.2018 | 0.7617 | 0.3878 | 0.8211 | 15.0120 | 1.3004 | |
| simple | 0.0508 | 0.2041 | 0.6952 | 0.3686 | 0.7183 | 14.9286 | 1.2711 | |
| t_exp | ideal | 0.0042 | 0.0083 | 0.0103 | 0.0099 | 0.0120 | 0.0231 | 0.0031 |
| mixed | 0.0358 | 0.1933 | 0.7128 | 0.3281 | 0.5861 | 11.2516 | 1.3474 | |
| simple | 0.0353 | 0.1931 | 0.6470 | 0.3241 | 0.5666 | 10.9228 | 1.2003 | |
| t_norm | ideal | 0.0037 | 0.0086 | 0.0107 | 0.0102 | 0.0126 | 0.0209 | 0.0031 |
| mixed | 10.6921 | 32.0753 | 58.9187 | 47.1701 | 67.8167 | 589.9308 | 53.0344 | |
| simple | 10.6071 | 28.9494 | 57.0330 | 46.4800 | 68.3398 | 593.1615 | 52.4424 | |
| t_snorm | ideal | 0.0042 | 0.0082 | 0.0102 | 0.0101 | 0.0121 | 0.0185 | 0.0029 |
| mixed | 0.1630 | 0.4268 | 0.7944 | 0.6313 | 0.9302 | 4.2049 | 0.5745 | |
| simple | 0.1562 | 0.4159 | 0.7559 | 0.6070 | 0.9392 | 3.0123 | 0.4905 |
Zwróćmy uwagę na średnią wartość MSE. Dla wszystkich przypadków model Ideal czyli lm(Y~X1+I(X2^2)) osiągał najniższe średnie wartości błędu MSE i najmniejszą zmienność błędu, co jest naturalne, ponieważ odpowiada on dokłądnie generującej funkcji. Błąd na poziomie 0.01 jest głównie spowodowany szumem losowym więc jest to najlepsze możliwe dopasowanie. W rzeczywistości jednak, raczej nie znamy funkcji generującej dane, a szukamy modelu, który najlepiej dopasuje się do danych. Tak traktujemy dwa pozostałe modele Mixed lm(Y~X1+X2+I(X1*X2)) i Simple lm(Y~X1+X2) . Dla tych modeli obserwujemy większe zróżnicowanie wartości MSE - między różnymi rozkładami, ale również w obrębie jednego przykładu. Przyjrzyjmy się kilku przypadkom, aby lepiej zrozumieć skąd te różnice się biorą.
# Wykres: średnie MSE (słupki)
ggplot(results_summary, aes(x = case, y = mean_MSE, fill = model)) +
geom_col(position = "dodge") +
theme_minimal() +
labs(title = "Średnie wartości MSE w poszczególnych przypadkach",
x = "Przypadek",
y = "Średnie MSE") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))Z wykresu widzimy, że dla 4 przypadków rozkładów modele mixed i simple poradziły sobie bardzo źle. Dla porównania zestawimy wyniki dla przypadków:
double_unif_narrow i double_unif_wide
t_norm i t_snorm
exp_norm i norm_exp
W tym celu wczytane zostaną dane z 1 ramki dla każdego przykładu
plot_case_fit <- function(case_name) {
df <- read.csv(paste0("data_samples/", case_name, ".csv"))
set.seed(123)
train_idx <- sample(1:nrow(df), size = 0.8 * nrow(df))
train <- df[train_idx, ]
test <- df[-train_idx, ]
m_simple <- lm(Y ~ X1 + X2, data = train)
m_ideal <- lm(Y ~ X1 + I(X2^2), data = train)
m_mixed <- lm(Y ~ X1 + X2 + I(X1 * X2), data = train)
test$pred_simple <- predict(m_simple, newdata = test)
test$pred_ideal <- predict(m_ideal, newdata = test)
test$pred_mixed <- predict(m_mixed, newdata = test)
mse_simple <- mean((test$Y - test$pred_simple)^2)
mse_ideal <- mean((test$Y - test$pred_ideal)^2)
mse_mixed <- mean((test$Y - test$pred_mixed)^2)
var_X1 <- var(test$X1)
var_X2 <- var(test$X2)
var_Y <- var(test$Y)
summary_table <- tibble(
case = case_name,
var_X1 = var_X1,
var_X2 = var_X2,
var_Y = var_Y,
mse_simple = mse_simple,
mse_ideal = mse_ideal,
mse_mixed = mse_mixed
)
# Wykresy
p1 <- ggplot(test, aes(x = X1, y = Y)) +
geom_point(alpha = 0.5, color = "grey40") +
geom_line(aes(y = pred_simple, color = "Simple"), linewidth = 0.8) +
geom_line(aes(y = pred_ideal, color = "Ideal"), linewidth = 0.8) +
geom_line(aes(y = pred_mixed, color = "Mixed"), linewidth = 0.8) +
labs(
title = paste("Dopasowanie względem X1"),
x = "X1", y = "Y", color = "Model"
) +
theme_minimal()
p2 <- ggplot(test, aes(x = X2, y = Y)) +
geom_point(alpha = 0.5, color = "grey40") +
geom_line(aes(y = pred_simple, color = "Simple"), linewidth = 0.8) +
geom_line(aes(y = pred_ideal, color = "Ideal"), linewidth = 0.8) +
geom_line(aes(y = pred_mixed, color = "Mixed"), linewidth = 0.8) +
labs(
title = paste("Dopasowanie względem X2"),
x = "X2", y = "Y", color = "Model"
) +
theme_minimal()
print(case_name)
grid.arrange(p1, p2, ncol = 2)
summary_table %>%
kbl(caption = paste("Statystyki dla przypadku:", case_name)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
}plot_pred_vs_actual <- function(case_name) {
df <- read.csv(paste0("data_samples/", case_name, ".csv"))
set.seed(123)
train_idx <- sample(1:nrow(df), size = 0.8 * nrow(df))
train <- df[train_idx, ]
test <- df[-train_idx, ]
m_simple <- lm(Y ~ X1 + X2, data = train)
m_ideal <- lm(Y ~ X1 + I(X2^2), data = train)
m_mixed <- lm(Y ~ X1 + X2 + I(X1 * X2), data = train)
# Predykcje
test$pred_simple <- predict(m_simple, newdata = test)
test$pred_ideal <- predict(m_ideal, newdata = test)
test$pred_mixed <- predict(m_mixed, newdata = test)
# Obliczanie MSE
mse_simple <- mean((test$Y - test$pred_simple)^2)
mse_ideal <- mean((test$Y - test$pred_ideal)^2)
mse_mixed <- mean((test$Y - test$pred_mixed)^2)
# Pomocnicza funkcja do tworzenia wykresu dla jednego modelu
make_plot <- function(pred_col, model_name, mse_value, color) {
ggplot(test, aes(x = Y, y = !!sym(pred_col))) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black") +
geom_segment(aes(xend = Y, yend = !!sym(pred_col)), alpha = 0.4, color = color) +
geom_point(color = color, size = 1.5) +
labs(
title = paste("Model", model_name),
subtitle = sprintf("MSE = %.4f", mse_value),
x = "Prawdziwy Y",
y = "Przewidziany Ŷ"
) +
theme_minimal() +
coord_equal()
}
# Trzy wykresy obok siebie
p1 <- make_plot("pred_simple", "Simple", mse_simple, "darkred")
p2 <- make_plot("pred_ideal", "Ideal", mse_ideal, "steelblue")
p3 <- make_plot("pred_mixed", "Mixed", mse_mixed, "darkgreen")
print(case_name)
grid.arrange(p1, p2, p3, nrow = 1)
}Obie wersje mają ten sam rodzaj rozkładu - jednostajny - ale różny zakres:
narrow: \(X_i \sim U(-3,3)\)
wide: \(X_i \sim U(-5,5)\)
plot_case_fit("double_unif_narrow")[1] "double_unif_narrow"
| case | var_X1 | var_X2 | var_Y | mse_simple | mse_ideal | mse_mixed |
|---|---|---|---|---|---|---|
| double_unif_narrow | 3.172489 | 3.013163 | 5.848907 | 3.165899 | 0.0081483 | 2.913976 |
plot_pred_vs_actual("double_unif_narrow")[1] "double_unif_narrow"
Na podstawie powyższego wykresu scatterplot możemy powiedzieć jak dobrze model dopasowuje się do danych. Im bliżej linii \(\hat{Y} = Y\) znajdują się punkty przedstawiające wartości przewidywane, tym lepszy jest model. Widać wyraźnie, że model Ideal “rzuca” przewidywania w większości pokrywające się z prawdziwymi wartościami.
diagnose_model <- function(case_name) {
df <- read.csv(paste0("data_samples/", case_name, ".csv"))
model <- lm(Y ~ X1 + I(X2^2), data = df)
par(mfrow = c(2, 2))
plot(model, main = paste("Diagnostyka modelu:", case_name))
par(mfrow = c(1, 1))
}diagnose_model("double_unif_narrow")Dla idealnie dopasowanego modelu wykres Residuals vs Fitted (reszty względem wartości dopasowanych) nie powinien wykazywać żadnego widocznego wzoru. Punkty powinny być rozmieszczone losowo wokół osi poziomej. Jeśli natomiast zauważylibyśmy np. kształt przypominający parabolę, mogłoby to oznaczać, że zależność między zmiennymi nie jest liniowa.
W naszym przypadku punkty są rozproszone nieregularnie, co sugeruje, że przyjęty model liniowy jest odpowiedni.
Wykres Normal Q-Q służy do oceny, czy reszty (błędy) mają rozkład normalny. Jeśli pochodzą z rozkładu normalnego, punkty powinny układać się wzdłuż prostej przerywanej.
W rozważanym przypadku obserwujemy, że punkty przebiegają wzdłuż tej linii, co potwierdza, że błędy mają rozkład zbliżony do normalnego.
Wykres Scale-Location (znany też jako Spread-Location) pokazuje, czy wariancja błędów jest stała (homoscedastyczność). Jeśli czerwona linia trendu jest prawie pozioma, a punkty są równomiernie rozproszone, możemy uznać, że założenie jednorodności wariancji jest spełnione.
Wykres Residuals vs Leverage pomaga zidentyfikować obserwacje wpływowe, czyli takie, które znacząco oddziałują na estymowane parametry modelu. Punkty znajdujące się poza liniami Cook’s distance (zaznaczonymi przerywanymi liniami) są potencjalnie wpływowe - ich usunięcie mogłoby istotnie zmienić dopasowanie modelu.
Nasz model nie posiada potencjalnie wpływowych obserwacji.
plot_case_fit("double_unif_wide")[1] "double_unif_wide"
| case | var_X1 | var_X2 | var_Y | mse_simple | mse_ideal | mse_mixed |
|---|---|---|---|---|---|---|
| double_unif_wide | 8.672501 | 9.898818 | 37.90714 | 24.04106 | 0.0102462 | 28.77621 |
plot_pred_vs_actual("double_unif_wide")[1] "double_unif_wide"
diagnose_model("double_unif_wide")Z powyższych wykresów możemy wyciągnąć podobne wnioski co w przypadku double_unif_narrow.
| “double_unif” model | narrow | wide |
|---|---|---|
| Simple | 3.121 | 23.947 |
| Ideal | 0.011 | 0.009 |
| Mixed | 2.858 | 28.41 |
Model idealny w obu przypadkach ma bardzo niskie MSE, a to oznacza, że poprawnie odwzorowuje funkcję generującą dane niezależnie od zakresu zmiennych.
Natomiast modele Simple i Mixed mają MSE znacznie większe dla „wide”, ponieważ wartości \(X_2\) i \(X_1\) są rozrzucone szerzej, co skutkuje większym zróżnicowaniem wartości \(Y\), a przez to większymi błędami predykcji - bo błędy kwadratowe \((Y - \hat{Y})^2\) rosną szybciej, im większy rozrzut (wariancja \(X_1\), \(X_2\)). Wartości \(Y\) też przyjmują większy zakres, bo są funkcją \(X_1\) i \(X_2\).
Wariancje zmiennych \(X_1\) i \(X_2\) w przypadkach narrow i wide różnią się kilkukrotnie. To oznacza, że skala błędów w modelach prostych (Simple, Mixed) rośnie proporcjonalnie do rozrzutu zmiennych, ponieważ MSE jest miarą bezwzględną, a nie znormalizowaną.
Sama skala i wariancja zmiennych wpływa na wartość MSE, nawet jeśli dopasowanie względne (np. \(R^2\)) jest podobne. Dlatego MSE nie jest bezwzględnie porównywalne między przypadkami o różnych skalach danych.
(Współczynnik determinacji \(R^2\) pokazuje, jak dobrze model wyjaśnia zmienność danych niezależnie od ich skali. MSE natomiast rośnie wraz ze skalą i wariancją zmiennych, dlatego przy porównywaniu modeli dla różnych rozkładów bardziej miarodajne jest porównanie \(R^2\) niż samych wartości MSE.
Porównujemy t_norm i t_snorm, dla których w obu przypadkach \(X_1 \sim t(3)\) czyli są to rozkłady t-Studenta z 3 stopniami swobody. Natomiast wariancja \(X_2\) się różni. t_snorm \(X_2 \sim N(0,1)\), a t_snorm \(X_2 \sim N(0,2.9)\)
Dla t_snorm obie zmienne mają umiarkowany rozrzut. Zmienna \(X_2\) w t_norm ma znacznie większe odchylenie standardowe, co oznacza, że dane są bardziej rozproszone w przestrzeni.
plot_case_fit("t_norm")[1] "t_norm"
| case | var_X1 | var_X2 | var_Y | mse_simple | mse_ideal | mse_mixed |
|---|---|---|---|---|---|---|
| t_norm | 2.310711 | 6.934366 | 20.48763 | 19.13618 | 0.0081343 | 23.31183 |
plot_pred_vs_actual("t_norm")[1] "t_norm"
diagnose_model("t_norm")plot_case_fit("t_snorm")[1] "t_snorm"
| case | var_X1 | var_X2 | var_Y | mse_simple | mse_ideal | mse_mixed |
|---|---|---|---|---|---|---|
| t_snorm | 1.183799 | 0.7243344 | 1.341954 | 0.2843283 | 0.0112532 | 0.3059577 |
Dla t_snorm wartości \(X_2\) są małe (blisko zera), więc składnik \(X_2^2\) nie ma dużego wpływu na \(Y\) -> Model łatwo się dopasowuje, błędy są małe.
Dla t_norm, wartości \(X_2\) są dużo większe - kwadrat rośnie szybko więc różnice w \(Y\) są dużo większe -> nawet mały błąd prognozy \(\hat{Y}\) przekłada się na duży błąd kwadratowy, dlatego MSE rośnie wykładniczo wraz z rozrzutem zmiennych.
| model | t_snorm | t_norm |
|---|---|---|
| Simple | 0.446 | 141.746 |
| Ideal | 0.01 | 0.01 |
| Mixed | 0.46 | 147.152 |
plot_pred_vs_actual("t_snorm")[1] "t_snorm"
diagnose_model("t_snorm")W t_norm zobaczymy większy rozrzut reszt, co wskazuje na heteroskedastyczność (rosnącą wariancję błędów przy większych wartościach X2).
Wariant t_norm ma znacznie większy rozrzut - większą wariancję \(X_2\) i cięższe ogony \(X_1\). Oznacza to, że w danych częściej pojawiają się wartości skrajne (outliery), przez co błędy kwadratowe przy modelach niedopasowanych (Simple, Mixed) są ekstremalnie duże.
Dla modelu Ideal różnice są niewielkie, bo model dokładnie odzwierciedla funkcję generującą dane. Dla modeli Simple i Mixed MSE dla t_norm jest kilkaset razy większe, mimo że dopasowanie względne (np. \(R^2\) ) może być podobne.
To dowód, że MSE jest wrażliwe na skalę danych i nie nadaje się do bezpośredniego porównywania między przypadkami o różnych wariancjach.
Wniosek: Rozkłady z większą wariancją lub cięższymi ogonami (np. t-Studenta, szeroki uniform) powodują znaczne zwiększenie MSE przy modelach uproszczonych, mimo identycznego modelu funkcji.
Przy norm_exp \(X_1 \sim N(0, 2.9)\) i \(X_2 \sim Exp(1.2)\). Więc jedna zmienna jest symetryczna, a druga skośna w prawo.
plot_case_fit("norm_exp")[1] "norm_exp"
| case | var_X1 | var_X2 | var_Y | mse_simple | mse_ideal | mse_mixed |
|---|---|---|---|---|---|---|
| norm_exp | 7.3134 | 1.201378 | 7.41053 | 0.2582772 | 0.0078207 | 0.2438117 |
plot_pred_vs_actual("norm_exp")[1] "norm_exp"
Oczywiście obserwujemy bardzo dobre dopasowanie modelu Ideal, ale co ciekawe modele Mixed i Simple mają małe MSE. Wynika to z zakresu i rozkładu \(X_2\):
\(X_2\) z rozkładu wykładniczego o \(\lambda =1.2\) przyjmuje małe wartości (średnio ok. 0.83), a większe wartości są bardzo rzadkie. W tym zakresie \(X_2^2 \approx X_2\) (bo przy małych wartościach, kwadrat nie zmienia dużo.
Czyli model \(Y \sim X_1 +X_2\) nie „widzi” mocno nieliniowości - i mimo że teoretycznie jest źle określony, praktycznie dopasowuje się bardzo dobrze, bo dane nie pokazują dużych odchyleń od liniowości.
Model \(Y \sim X_1 +X_2+X1\cdot X2\) osiągnął wartość MSE bardzo podobną do Simple. Składnik \(X1\cdot X2\) jest niepotrzebny, jednak nie ma on dużego wpływu na model, jeśli jego współczynnik jest mały.
diagnose_model("norm_exp")W przypadku exp_norm \(X_1 \sim Exp(1.2)\) i \(X_2 \sim N(0, 2.9)\)
plot_case_fit("exp_norm")[1] "exp_norm"
| case | var_X1 | var_X2 | var_Y | mse_simple | mse_ideal | mse_mixed |
|---|---|---|---|---|---|---|
| exp_norm | 1.114356 | 9.109857 | 107.2702 | 108.5138 | 0.0079043 | 124.1917 |
plot_pred_vs_actual("exp_norm")[1] "exp_norm"
Gdy zamienimy rozkłady, z których pochodzą zmienne objaśniające, drastycznie zwiększa się wartość błędu średniokwadratowego dla modeli Simple i Mixed.
W tym wypadku składnik \(0.6\cdot X_2^2\) ma szeroki rozrzut. Składnik \(0.6\cdot X_2^2\) dominuje w Y, czego prosty model liniowy nie jest w stanie uchwycić.
diagnose_model("exp_norm")Model nie w pełni dopasował się do liniowej zależności w danych - punkty na 1. wykresie są skupione w okolicach 0.
Dla dużych wartości przewidywanych model ma większy rozrzut błędów -> heteroscedastyczność. Model lepiej radzi sobie dla małych Y, gorzej dla dużych.
Większość punktów ma bardzo małą dźwignię, model nie posiada obserwacji wpływowych.
Oznacza to, że zamiana rozkładów zmiennych objaśniających może diametralnie zmienić charakter zależności i dopasowanie modeli, nawet jeśli sama struktura wzoru pozostaje taka sama.
Pani Profesor Julie Josse na zimowej szkole z Causality i XAI w Srobonne University, przedstawiła problem zmiennych ukrytych w bardzo intuicyjny i zabawny sposób aby zrozumieć dlaczego jest to istotne i czemu same mlowe i statystyczne podejście może nas zawieść. Możemy np. dowieść statystycznie, że spanie w butach powoduje ból głowy następnego dnia. Pani Profesor dowiodła to testem statystycznym z danych, że ta korelacja jest istotna. W danych nie mieliśmy jednak informacji o tym, że dana osoba spożyła wcześniej spore ilości alkoholu, który sprawił, że po powrocie do domu z imprezy nie była w stanie zdjąć butów i poszła spać. Rano obudziła się z kacem - tego nam same dane nie powiedzą i potrzebujemy causal approach. Posłużyłem się chatem gpt do wygenerowania poniższej grafiki w celu ilustracji, abyśmy zrozumieli lepiej graf przyczynowo skutkowo korelacyjny (chat nie mógł zrozumieć prawidłowo kierunku strzałek więc zostawiam grafikę tak jak ja wygenerował :D).
Przechodząc do konkretnego przypadku naszego projektu w tym podrozdziale badamy przypadek, w którym nieobserwowana zmienna \(Z\) oddziałuje zarówno na zmienne \(X1, X2\)jak i na zmienną objaśnianą \(Y\) Taki układ powoduje:
zależność pomiędzy \(X1\) i \(X2\)
systematyczny składnik wpływający na \(Y\) którego nie potrafimy „wyjaśnić” przez posiadane zmienne z macierzy eksperymentu, tzw. bias przez pominięcie
potencjalnie złudną poprawę na danych uczących i gorszą generalizację.
Ponieważ \(Z\) nie jest w rozpatrywanej ramce danych, pokażemy, jak jego wpływ widoczny jest pośrednio: w strukturze predyktorów oraz w resztach modeli oraz w zmianie MSE.
Zbadajmy to jaki “cień” zostawiła w naszych zmiennych \(X-owych\) nieznana zmienna \(Z\) za pomocą techniki PCA (analizy głównych składowych) a jak wiemy główną składową każdej z zmiennych \(X-owych\) jest ukryta i bardzo tajemnicza zmienna \(Z\)
Czy widać tak naprawdę wpływ zmiennej \(Z\) na zmienne \(X-owe\)?
library(dplyr)
library(tibble)
library(purrr)
library(ggplot2)
library(knitr)
library(kableExtra)
library(scales)
dfs <- list(
df1 = df1_hidden_impact,
df2 = df2_hidden_impact,
df3 = df3_hidden_impact
)pc1_score <- function(data){ pr <- prcomp(scale(dplyr::select(data, X1, X2)))
list(score = as.numeric(pr$x[,1]), var_expl = pr$sdev[1]^2 / sum(pr$sdev^2)) }
proxy_tbl <- imap_dfr(dfs, function(dfi, nm){
pc1 <- pc1_score(dfi)
tibble(
zbior = nm,
cor_X1_X2 = cor(dfi$X1, dfi$X2),
var_expl_PC1 = round(pc1$var_expl, 3),
cor_PC1_X1 = cor(pc1$score, dfi$X1),
cor_PC1_X2 = cor(pc1$score, dfi$X2)
)
})
proxy_tbl %>%
transmute(
`Zbiór` = zbior,
`Corr(X1, X2)` = round(cor_X1_X2, 3),
`Corr(PC1, X1)` = round(cor_PC1_X1, 3),
`Corr(PC1, X2)` = round(cor_PC1_X2, 3),
`PC1: udział wariancji`= scales::percent(var_expl_PC1, accuracy = 0.1)
) %>%
kable(align = "lcccc",
booktabs = TRUE,
caption = "Sygnały wspólnej przyczyny (proxy-Z)") %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover")) %>%
column_spec(1, bold = TRUE) %>%
add_header_above(c(" " = 1, "Korelacje" = 3, " " = 1))| Zbiór | Corr(X1, X2) | Corr(PC1, X1) | Corr(PC1, X2) | PC1: udział wariancji |
|---|---|---|---|---|
| df1 | 1 | -1 | -1 | 100.0% |
| df2 | 1 | 1 | 1 | 100.0% |
| df3 | 1 | -1 | -1 | 100.0% |
Oczywiście że tak ponieważ w wariancie hidden_impact przyjęliśmy deterministyczne zależności. Skutkiem jest rząd macierzy cech równy 1, tj. \(X1\) i \(X2\) leżą na jednej prostej, co daje bardzo silne zależności, PC1 odtwarza kierunek \(Z\) z dokładnością do znaku. Jest to zamierzony, skrajny przykład confoundingu: cała zmienność predyktorów pochodzi ze wspólnej, nieobserwowanej przyczyny.
pc1_all <- purrr::imap_dfr(dfs, function(dfi, nm){
pc <- pc1_score(dfi)
tibble(
zbior = nm,
PC1 = pc$score,
var_expl = pc$var_expl
)
})
facet_labels <- pc1_all %>%
distinct(zbior, var_expl) %>%
mutate(label = paste0(zbior, " - PC1: ", percent(var_expl, accuracy = 0.1))) %>%
{ setNames(.$label, .$zbior) }
pc1_hist <- ggplot(pc1_all, aes(x = PC1)) +
geom_histogram(aes(y = ..density..), bins = 30, alpha = 0.75) +
geom_density(adjust = 1.2, linewidth = 0.8) +
geom_vline(xintercept = 0, linetype = "dashed") +
facet_wrap(~ zbior, scales = "free_x", labeller = as_labeller(facet_labels)) +
labs(
title = "Rozkład PC1 (proxy-Z) w zestawach hidden_impact",
x = "PC1 (standaryzowany score z PCA na {X1, X2})",
y = "Gęstość"
) +
theme_minimal(base_size = 12) +
theme(
strip.text = element_text(face = "bold"),
plot.title = element_text(face = "bold")
)
pc1_histWniosek z wykresu jest widoczny gołym okiem ponieważ 100% udziału wskazuje że cała zmienność zmiennych \(X1\) i \(X2\) leży w jednym kierunku. Co potwierdza wniosek z poprzedniej powyższej tabeli.
Gdy PC1 wyjaśnia bardzo dużo wariancji, mamy silne przesłanki, że w danych działa nieobserwowany czynnik.
W takiej sytuacji modele regresji bez tego czynnika będą zostawiać strukturę w resztach (co pokazujemy w kolejnych wykresach) i będą miały podwyższony błąd generalizacji.
To uzasadnia, dlaczego w praktyce warto szukać proxy (zmiennych zastępczych) lub instrumentów żeby zredukować wpływ pominiętej zmiennej.
Ustalmy dwie formy modelu którymi będziemy się posługiwać w dalszych rozważaniach
\[ Y \sim X1 + X2 \]
\[ Y \sim X1 + I(X_2^2) \]
form_simple <- Y ~ X1 + X2
form_ideal <- Y ~ X1 + I(X2^2)cv_mse <- function(form, data, K = 5, seed = 123){
set.seed(seed)
n <- nrow(data)
folds <- sample(rep(1:K, length.out = n)) # losowe przypisanie do K foldów
mse_folds <- sapply(1:K, function(k){
tr <- data[folds != k, , drop = FALSE]
te <- data[folds == k, , drop = FALSE]
fit <- lm(form, data = tr)
mean((te$Y - predict(fit, newdata = te))^2)
})
list(mean = mean(mse_folds),
sd = sd(mse_folds),
se = sd(mse_folds) / sqrt(K),
per_fold = mse_folds)
}K <- 5
seed_cv <- 123
cv_tbl <- imap_dfr(dfs, function(dfi, nm){
cv_s <- cv_mse(form_simple, dfi, K = K, seed = seed_cv)
cv_i <- cv_mse(form_ideal, dfi, K = K, seed = seed_cv)
tibble(
zbior = nm,
model = c("simple", "ideal"),
CV_MSE = c(cv_s$mean, cv_i$mean),
SE = c(cv_s$se, cv_i$se)
)
})
cv_folds_long <- imap_dfr(dfs, function(dfi, nm){
list(
tibble(zbior = nm, model = "simple", fold = 1:K,
mse = cv_mse(form_simple, dfi, K = K, seed = seed_cv)$per_fold),
tibble(zbior = nm, model = "ideal", fold = 1:K,
mse = cv_mse(form_ideal, dfi, K = K, seed = seed_cv)$per_fold)
) |> bind_rows()
})library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(scales)
table_pretty <- cv_tbl %>%
dplyr::select(zbior, model, CV_MSE, SE) %>%
dplyr::mutate(CI_low = CV_MSE - 1.96 * SE,
CI_hi = CV_MSE + 1.96 * SE) %>%
tidyr::pivot_wider(names_from = model, values_from = c(CV_MSE, SE, CI_low, CI_hi)) %>%
dplyr::mutate(
`Δ (simple - ideal)` = CV_MSE_simple - CV_MSE_ideal,
`Rel. poprawa ideal vs simple` = ifelse(
is.finite(CV_MSE_simple) & CV_MSE_simple != 0,
(CV_MSE_simple - CV_MSE_ideal) / CV_MSE_simple, NA_real_
)
)
table_pretty %>%
dplyr::mutate(
dplyr::across(
.cols = c(CV_MSE_simple, SE_simple, CI_low_simple, CI_hi_simple,
CV_MSE_ideal, SE_ideal, CI_low_ideal, CI_hi_ideal,
`Δ (simple - ideal)`),
~ round(.x, 4)
),
`Rel. poprawa ideal vs simple` = percent(`Rel. poprawa ideal vs simple`, accuracy = 0.1)
) %>%
dplyr::rename(
Zbiór = zbior,
`CV-MSE` = CV_MSE_simple, `SE` = SE_simple,
`CI (low)` = CI_low_simple, `CI (hi)` = CI_hi_simple,
`CV-MSE ` = CV_MSE_ideal, `SE ` = SE_ideal,
`CI (low) ` = CI_low_ideal, `CI (hi) ` = CI_hi_ideal
) %>%
knitr::kable(align = "lrrrrrrrrrr", booktabs = TRUE,
caption = "K-fold CV-MSE (K=5): porównanie modeli bez Z (średnia, SE, 95% CI, różnice)") %>%
kableExtra::kable_styling(full_width = FALSE, bootstrap_options = c("striped","hover")) %>%
kableExtra::add_header_above(c(" " = 1, "Model: simple" = 4, "Model: ideal" = 4, "Różnice" = 2)) %>%
kableExtra::column_spec(1, bold = TRUE)| Zbiór | CV-MSE | CV-MSE | SE | SE | CI (low) | CI (low) | CI (hi) | CI (hi) | Δ (simple - ideal) | Rel. poprawa ideal vs simple |
|---|---|---|---|---|---|---|---|---|---|---|
| df1 | 4.5094 | 0.0126 | 1.4715 | 0.0020 | 1.6253 | 0.0086 | 7.3935 | 0.0165 | 4.4968 | 99.7% |
| df2 | 2.1652 | 0.0095 | 0.5558 | 0.0009 | 1.0759 | 0.0077 | 3.2545 | 0.0113 | 2.1557 | 99.6% |
| df3 | 2.5123 | 0.0105 | 1.0422 | 0.0016 | 0.4695 | 0.0074 | 4.5551 | 0.0135 | 2.5018 | 99.6% |
Gdzie:
CV_MSE_model - średni błąd walidacji krzyżowej (K-fold) danego modelu.
SE_model – błąd standardowy średniej CV-MSE (zmienność między foldami).
CI_low_model, CI_hi_model – dolna i górna granica 95% CI dla CV-MSE (≈ CV_MSE ± 1.96·SE).
Specyfikacja modelu ma pierwszorzędne znaczenie. Sama obecność nieobserwowanej \(Z\) nie musi niszczyć jakości predykcji, jeśli informacja o \(Z\) jest „pośrednio zakodowana” w \(X\) i uchwycona właściwą postacią
Brak właściwej postaci ma wpływ duży koszt generalizacji. simple traci nie dlatego, że nie widzi \(Z\), tylko dlatego, że źle modeluje zależność \(Y\) do \(X\).
Konkluzja. W naszym skrajnym scenariuszu hidden_impact dominująca struktura w \(X\) (PC1 ≈ cień \(Z\)) oraz CV-MSE pokazują, że kluczowa jest prawidłowa postać funkcjonalna względem obserwowalnych cech; dopiero gdy po jej uwzględnieniu reszty nadal są skorelowane z proxy-Z i CV-MSE pozostaje wysokie, sensowne jest inwestowanie w dodatkowe zmienne/proxy/instrumenty.
plot_case_fit2 <- function(case_name) {
df <- read.csv(paste0("data_samples/", case_name, ".csv"))
set.seed(123)
train_idx <- sample(1:nrow(df), size = 0.8 * nrow(df))
train <- df[train_idx, ]
test <- df[-train_idx, ]
m_simple <- lm(Y ~ X1 + X2, data = train)
m_ideal <- lm(Y ~ X1 + I(X2^2), data = train)
m_mixed <- lm(Y ~ X1 + X2 + I(X1 * X2), data = train)
test$pred_simple <- predict(m_simple, newdata = test)
test$pred_ideal <- predict(m_ideal, newdata = test)
test$pred_mixed <- predict(m_mixed, newdata = test)
mse_simple <- mean((test$Y - test$pred_simple)^2)
mse_ideal <- mean((test$Y - test$pred_ideal)^2)
mse_mixed <- mean((test$Y - test$pred_mixed)^2)
var_X1 <- var(test$X1)
var_X2 <- var(test$X2)
var_Y <- var(test$Y)
summary_table <- tibble(
case = case_name,
var_X1 = var_X1,
var_X2 = var_X2,
var_Y = var_Y,
mse_simple = mse_simple,
mse_ideal = mse_ideal,
mse_mixed = mse_mixed
)
# Wykresy
p1 <- ggplot(test, aes(x = X1, y = Y)) +
geom_point(alpha = 0.5, color = "grey40") +
geom_line(aes(y = pred_simple, color = "Simple"), linewidth = 0.8) +
geom_line(aes(y = pred_ideal, color = "Ideal"), linewidth = 0.8) +
#geom_line(aes(y = pred_mixed, color = "Mixed"), linewidth = 0.8) +
labs(
title = paste("Dopasowanie względem X1"),
x = "X1", y = "Y", color = "Model"
) +
theme_minimal()
p2 <- ggplot(test, aes(x = X2, y = Y)) +
geom_point(alpha = 0.5, color = "grey40") +
geom_line(aes(y = pred_simple, color = "Simple"), linewidth = 0.8) +
geom_line(aes(y = pred_ideal, color = "Ideal"), linewidth = 0.8) +
#geom_line(aes(y = pred_mixed, color = "Mixed"), linewidth = 0.8) +
labs(
title = paste("Dopasowanie względem X2"),
x = "X2", y = "Y", color = "Model"
) +
theme_minimal()
print(case_name)
grid.arrange(p1, p2, ncol = 2)
summary_table %>%
kbl(caption = paste("Statystyki dla przypadku:", case_name)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
}
plot_case_fit2("hidden_impact")[1] "hidden_impact"
| case | var_X1 | var_X2 | var_Y | mse_simple | mse_ideal | mse_mixed |
|---|---|---|---|---|---|---|
| hidden_impact | 9.894365 | 2.660213 | 1118.153 | 5.524909 | 0.0136601 | 0.0136601 |
diagnose_model("hidden_impact")Zdefiniowano “prawdziwą” funkcję generującą Y oraz globalny szum (ten sam poziom dla wszystkich scenariuszy), co umożliwia porównywanie MSE między przypadkami.
Rozkłady predyktorów: normalne (różne wariancje), jednostajne (wąskie/szerokie), wykładnicze, Poissona, t-Studenta, mieszane (w tym bimodalne) i kombinacje asymetryczne/ciężkoogonowe.
Budowano modele dopasowane do struktury danych (proste vs „idealne” względem mechanizmu) i porównywano CV-MSE (K=5) wraz z przedziałami ufności.
Globalny szum utrzymywał wspólną skalę błędu, więc różnice w MSE wynikały głównie z rozkładów, postaci funkcjonalnej i obecności/siły Z.
Rozkłady skośne i ciężkoogonowe (t-Student, wykładniczy, mieszane) zwiększają MSE i jego zmienność między foldami; modele są bardziej czułe na obserwacje skrajne.
Niedopasowanie postaci funkcjonalnej (np. pominięcie nieliniowości) podnosi MSE silniej niż sama zmiana rozkładu – to koszt uogólniania.
Standaryzacja i rozsądny wybór featurów stabilizują MSE, ale nie kompensują błędnej postaci modelu.
Z wprowadzono jako wspólną przyczynę, konstruowano proxy-Z (np. komponent główny).
Specyfikacja modelu decyduje – poprawna postać potrafi „pośrednio” uchwycić wpływ Z i obniżyć MSE nawet bez jawnego Z.
Gdy wpływ Z jest dominujący, warto inwestować w dodatkowe zmienne/proxy/instrumenty; gdy dominuje błąd postaci, najpierw naprawiamy formę modelu.
Różnice MSE między modelem prostym a „idealnym” są czytelne i stabilne, co potwierdza powyższe tezy.
Diana Morzhak – pokazała, że wrażliwość MSE na strukturę danych jest użyteczną analogią do mechanizmów predictive coding
Mateusz Walo – przeprowadził causal approach: izolowanie efektu nieobserwowalnego czynnika, budowa proxy-Z, nadzorowanie i koordynowania prac i terminów realizacji projektu
Dominika Zydorczyk – usystematyzowała interpretację i komunikację wyników MSE dla różnych rozkładów/błędów, pokazując jak metryka odzwierciedla dopasowanie i kiedy myli modelarza.
Zespół w krótkim czasie wykonał dużo pracy: od przemyślanego generatora danych, przez rzetelną walidację, po klarowne wnioski – każdy dowiózł swój epsilon progres, a razem osiągnęliśmy dojrzałe rozumienie, jak i dlaczego MSE zachowuje się tak, jak obserwowaliśmy.